library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.0
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(shiny)
library(shinydashboard)
##
## Attaching package: 'shinydashboard'
##
## The following object is masked from 'package:graphics':
##
## box
library(ggplot2)
library(naniar)
library(lubridate)
## Loading required package: timechange
##
## Attaching package: 'lubridate'
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(dplyr)
library(RColorBrewer)
library(paletteer)
library(readr)
library(ggmap)
## ℹ Google's Terms of Service: <]8;;https://mapsplatform.google.comhttps://mapsplatform.google.com]8;;>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(devtools)
## Loading required package: usethis
library(here)
## here() starts at C:/GitHub Repositories/GOBY - Edits
library(ggrepel)
visibility <- read_csv("data/CA_visibility_data.csv")
## Rows: 26633 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (11): Dataset, SiteCode, Date, SiteName, State, ammNO3f_Unit, ammSO4f_Un...
## dbl (11): POC, Percentile, Latitude, Longitude, Elevation, ammNO3f_Val, ammS...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
visibility_tidy <- visibility%>%
na_if(-999)
Step 1: Create a Base Map
visibility_tidy %>%
select(Latitude, Longitude) %>% #note that the lat and long are in degrees decimal
summary()
## Latitude Longitude
## Min. :33.46 Min. :-124.1
## 1st Qu.:34.73 1st Qu.:-121.2
## Median :37.22 Median :-119.7
## Mean :37.32 Mean :-119.7
## 3rd Qu.:38.98 3rd Qu.:-118.1
## Max. :41.71 Max. :-116.4
lat <- c(33.46, 41.71)
long <- c(-124.1, -116.4)
bbox <- make_bbox(long, lat, f = 0.05)
map1 <- get_map(bbox, maptype = "terrain", source = "stamen")
## ℹ Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.
ggmap(map1)
Step 2: Adding Points to Base Map
I. Basic map with site locations
ggmap(map1) +
geom_point(data = visibility_tidy, aes(Longitude, Latitude), size=1.5) +
labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")
#ggsave("plot.png", width = 5, height = 5)
II.Basic map with site location and by elevation
ggmap(map1) +
geom_point(data = visibility_tidy, aes(Longitude, Latitude, colour=Elevation), size=2.5) +
labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")
III.Basic map with site location and colour name key - not informative
ggmap(map1) +
geom_point(data = visibility_tidy, aes(Longitude, Latitude, colour=SiteName), size=2.5) +
labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")
IV.Basic map with only ONE site location - not informative
visibility_yos <- visibility_tidy %>%
filter(SiteName=="Yosemite NP")
ggmap(map1) +
geom_point(data = visibility_yos, aes(Longitude, Latitude), size=2.5) +
labs(x= "Longitude", y= "Latitude", title="Yosemite Air Quality Research Site")
Saving maps
ggsave("plot.png", width = 5, height = 5) #Saves last plot as 5’ x 5’ file named "plot.png" in working directory.
#Matches file type to file extension.
Split date by year:
visibility2 <- visibility_tidy %>%
separate(Date, into= c("month", "day", "year"), sep = "/")
visibility3 <- visibility2 %>%
group_by(SiteName, year) %>%
summarise(across(c(SVR_Val, ammNO3f_Val, ammSO4f_Val, ECf_Val, OMCf_Val, SOILf_Val, Longitude, Latitude), mean, na.rm=T), total=n())
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
A) Graph Mean Standard Visual Range over time, faceted by site
visibility3 %>%
ggplot(aes(x=year, y=SVR_Val))+
geom_point(na.rm=TRUE)+
facet_wrap(~SiteName, ncol=4)+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
labs(title = "Standard Visual Range over Time",
x = "Year",
y = "Average Standard Visual Range (km)",) #fill by season?
B) Graph Mean Ammonium Nitrate Concentration Over Time faceted by site
visibility3 %>%
ggplot(aes(x=year, y=ammNO3f_Val))+
geom_point(na.rm=TRUE)+
facet_wrap(~SiteName, ncol=4)+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
labs(title = "Ammonium Nitrate Concentration over Time",
x = "Year",
y = "Ammonium Nitrate Concentration (ug/m^3)",) #fill by season?
C) Graph Mean Ammonium Sulfate Concentration Over Time faceted by site
visibility3 %>%
ggplot(aes(x=year, y=ammSO4f_Val))+
geom_point(na.rm=TRUE)+
facet_wrap(~SiteName, ncol=4)+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
labs(title = "Ammonium Sulfate Concentration over Time",
x = "Year",
y = "Ammonium Sulfate Concentration (ug/m^3)",) #fill by season?
D) Graph Mean Elemental Carbon Concentration Over Time faceted by site
visibility3 %>%
ggplot(aes(x=year, y=ECf_Val))+
geom_point(na.rm=TRUE)+
facet_wrap(~SiteName, ncol=4)+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
labs(title = "Elemental Carbon Concentration over Time",
x = "Year",
y = "Elemental Carbon Concentration (ug/m^3)",) #fill by season?
Notes: Interesting fluctuations for Lassen, Lava Beds, Point Reyes, Yosemite, Sequoia, Trinity. Link to wildfire data?
Separate Dates for plotting
visibility_mdy <- visibility_tidy%>%
separate(Date, into = c("month", "day", "year"))
visibility_mdy
## # A tibble: 26,633 × 24
## Dataset SiteCode POC month day year Percentile SiteName Latit…¹ Longi…²
## <chr> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl>
## 1 IMPRHR2 AGTI1 1 1 3 2011 10 Agua Tib… 33.5 -117.
## 2 IMPRHR2 AGTI1 1 1 6 2011 10 Agua Tib… 33.5 -117.
## 3 IMPRHR2 AGTI1 1 1 9 2011 70 Agua Tib… 33.5 -117.
## 4 IMPRHR2 AGTI1 1 1 12 2011 10 Agua Tib… 33.5 -117.
## 5 IMPRHR2 AGTI1 1 1 15 2011 10 Agua Tib… 33.5 -117.
## 6 IMPRHR2 AGTI1 1 1 18 2011 10 Agua Tib… 33.5 -117.
## 7 IMPRHR2 AGTI1 1 1 21 2011 10 Agua Tib… 33.5 -117.
## 8 IMPRHR2 AGTI1 1 1 24 2011 10 Agua Tib… 33.5 -117.
## 9 IMPRHR2 AGTI1 1 1 27 2011 10 Agua Tib… 33.5 -117.
## 10 IMPRHR2 AGTI1 1 1 30 2011 50 Agua Tib… 33.5 -117.
## # … with 26,623 more rows, 14 more variables: Elevation <dbl>, State <chr>,
## # ammNO3f_Val <dbl>, ammNO3f_Unit <chr>, ammSO4f_Val <dbl>,
## # ammSO4f_Unit <chr>, ECf_Val <dbl>, ECf_Unit <chr>, OMCf_Val <dbl>,
## # OMCf_Unit <chr>, SOILf_Val <dbl>, SOILf_Unit <chr>, SVR_Val <dbl>,
## # SVR_Unit <chr>, and abbreviated variable names ¹Latitude, ²Longitude
Make Plots Yosemite National Park - 2016 -> 2017 -> 2018 Empire Fire, Ferguson Fire
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Yosemite NP") %>%
summarize(mean_ano3_year= mean(ammNO3f_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_ano3_year)) +
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Ferguson Fire",
x = 8,
y = .59,
color = "black",
fill = NA) +
labs(title = "Mean Ammonium Nitrate in Yosemite National Park",
x = "Year",
y = "Ammonium Nitrate (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("yosemitenp_ammno3.png",
plot = last_plot())
## Saving 7 x 5 in image
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Yosemite NP") %>%
summarize(mean_aso4_year= mean(ammSO4f_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_aso4_year))+
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Ferguson Fire",
x = 8,
y = .79,
color = "black",
fill = NA) +
labs(title = "Mean Ammonium Sulfate in Yosemite National Park",
x = "Year",
y = "Ammonium Sulfate (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("yosemitenp_ammso4.png",
plot = last_plot())
## Saving 7 x 5 in image
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Yosemite NP") %>%
summarize(mean_omc_year= mean(OMCf_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_omc_year))+
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Ferguson Fire",
x = 8,
y = 10,
color = "black",
fill = NA) +
labs(title = "Mean OMC in Yosemite National Park",
x = "Year",
y = "Organic Mass by Carbon (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("yosemitenp_omc.png",
plot = last_plot())
## Saving 7 x 5 in image
Lassen Volcanic National Park - 2020 -> 2021 Dixie Fire
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Lassen Volcanic NP") %>%
summarize(mean_ano3_year= mean(ammNO3f_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_ano3_year)) +
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Dixie Fire",
x = 11,
y = .22,
color = "black",
fill = NA) +
labs(title = "Mean Ammonium Nitrate in Lassen Volcanic National Park",
x = "Year",
y = "Ammonium Nitrate (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("lassenvolcanicnp_ammno3.png",
plot = last_plot(),
width = 7,
height = 5)
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Lassen Volcanic NP") %>%
summarize(mean_aso4_year= mean(ammSO4f_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_aso4_year))+
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Dixie Fire",
x = 11,
y = .465,
color = "black",
fill = NA) +
labs(title = "Mean Ammonium Sulfate in Lassen Volcanic National Park",
x = "Year",
y = "Ammonium Sulfate (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("lassenvolcanicnp_ammso4.png",
plot = last_plot(),
width = 7,
height = 5)
visibility_mdy %>%
group_by(SiteName, year) %>%
filter(SiteName=="Lassen Volcanic NP") %>%
summarize(mean_omc_year= mean(OMCf_Val, na.rm = T)) %>%
ggplot(aes(x=year, y=mean_omc_year))+
geom_line(color="red",
alpha = 0.8,
group = 1) +
geom_point() +
geom_label(label = "Dixie Fire",
x = 11,
y = 5.6,
color = "black",
fill = NA) +
labs(title = "Mean OMC in Lassen Volcanic National Park",
x = "Year",
y = "Organic Mass by Carbon (ug/m^3)")
## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
ggsave("lassenvolcanicnp_omc.png",
plot = last_plot(),
width = 7,
height = 5)
improve_visibility <- visibility %>%
clean_names() %>%
na_if("-999") %>%
select(date, site_name, state,svr_val) %>%
separate(date, into = c("month", "day", "year"))
ui <- dashboardPage(
dashboardHeader(title = tags$div(style = "font-size: 30px;",
"Visual Range CA")),
dashboardSidebar(disable = TRUE),
dashboardBody(
tags$head(
tags$style(HTML(".box{border-radius: 20px; box-shadow: 5px 5px 10px #888888;}"))),
fluidRow( class="text center",
box(title= "Plot Options",width = 2,
radioButtons("x", "Choose The Site", choices = c("Agua Tibia", "Bliss SP (TRPA)", "Bliss SP (TRPA) (RHTS)", "Death Valley NP", "Dome Lands Wilderness", "Fresno", "Hoover", "Joshua Tree NP", "Kaiser", "Lake Tahoe Community College", "Lassen Volcanic NP", "Lava Beds NM", "Owens Valley", "Pinnacles NM", "Point Reyes National Seashore", "Redwood NP", "San Gabriel", "San Gorgonio Wilderness", "San Rafael", "Sequoia NP", "Trinity", "Wrightwood", "Yosemite NP")
)
), #Close box
box(title = h2("Standard Visual Range Over the Year", align="center"), width = 6, align="center",solidHeader = T,
plotOutput("plot", width = "550px", height = "525px")),
box(title="", width = 4, height = 650,
imageOutput("image"), style ="text-align: center;"
),
) #close fluidrow
) #close dashboardbody
) #close dashboard page
server <- function(input, output, session) {
session$onSessionEnded(stopApp)
output$plot <- renderPlot({
improve_visibility %>%
group_by(site_name, year) %>%
summarize("Mean_Standard_Visual_Range_km"= mean(svr_val, na.rm = T)) %>%
filter(site_name==input$x) %>%
ggplot(aes(x=year, y=Mean_Standard_Visual_Range_km,group=2, color=Mean_Standard_Visual_Range_km))+
geom_line(size=1.5,alpha=0.8)+
geom_point(size=5)+
theme(axis.title.x = element_text(size = 23, margin = margin(t = 15)))+
theme(axis.title.y = element_text(size = 23, margin = margin(r = 20)))+
theme(axis.text = element_text(size = 14))+
labs(y = "Average Standard Visual Range (km)")+
labs(x = "Year")+
scale_color_gradient(low="blue", high="red")+
guides(color = FALSE)
})
output$image <- renderImage({
list(src = "map.png", height = "570px", width="450px")
}, deleteFile = FALSE)
}
shinyApp(ui, server, options = list(launch.browser = TRUE))
Designating Seasons Warning: This code block is not time efficient (O(n)), expect a long run time.
seasons <- c()
for (date in visibility_tidy$Date){
date <- mdy(date)
# Winter Dec 20 - Mar 19
if ((month(date) == 12 & day(date) >= 20) | month(date) %in% c(1, 2) | (month(date) == 3 & day(date) <= 19)){
seasons <- c(seasons, "winter")
}
# Spring Mar 20 - June 20
else if ((month(date) == 3 & day(date) >= 20) | month(date) %in% c(4, 5) | (month(date) == 6 & day(date) <= 20)){
seasons <- c(seasons, "spring")
}
# Summer June 21 - Sept 22
else if ((month(date) == 6 & day(date) >= 21) | month(date) %in% c(7, 8) | (month(date) == 9 & day(date) <= 22)){
seasons <- c(seasons, "summer")
}
# Autumn Sept 23 - Dec 19
else{
seasons <- c(seasons, "autumn")
}
}
visibility_tidy$season = seasons
Separate Dates for plotting
visibility_tidy_mdy <- visibility_tidy%>%
separate(Date, into = c("month", "day", "year"))
head(visibility_tidy_mdy)
## # A tibble: 6 × 25
## Dataset SiteCode POC month day year Percentile SiteName Latit…¹ Longi…²
## <chr> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl>
## 1 IMPRHR2 AGTI1 1 1 3 2011 10 Agua Tibia 33.5 -117.
## 2 IMPRHR2 AGTI1 1 1 6 2011 10 Agua Tibia 33.5 -117.
## 3 IMPRHR2 AGTI1 1 1 9 2011 70 Agua Tibia 33.5 -117.
## 4 IMPRHR2 AGTI1 1 1 12 2011 10 Agua Tibia 33.5 -117.
## 5 IMPRHR2 AGTI1 1 1 15 2011 10 Agua Tibia 33.5 -117.
## 6 IMPRHR2 AGTI1 1 1 18 2011 10 Agua Tibia 33.5 -117.
## # … with 15 more variables: Elevation <dbl>, State <chr>, ammNO3f_Val <dbl>,
## # ammNO3f_Unit <chr>, ammSO4f_Val <dbl>, ammSO4f_Unit <chr>, ECf_Val <dbl>,
## # ECf_Unit <chr>, OMCf_Val <dbl>, OMCf_Unit <chr>, SOILf_Val <dbl>,
## # SOILf_Unit <chr>, SVR_Val <dbl>, SVR_Unit <chr>, season <chr>, and
## # abbreviated variable names ¹Latitude, ²Longitude
Explore the functionality of separating by season Finding the mean data per year per season.
mean_by_season <- visibility_tidy_mdy%>%
mutate(SiteName=paste(SiteName, year, sep="_"))%>%
mutate(SiteName=paste(SiteName, season, sep=" "))%>%
group_by(SiteName)%>%
summarise("mean_ammNO3_ppb"=mean(ammNO3f_Val, na.rm=T),
"mean_ammSO4_ppb"=mean(ammSO4f_Val, na.rm=T),
"mean_EC_ppb"=mean(ECf_Val, na.rm=T),
"mean_OMC_ppb"=mean(OMCf_Val, na.rm=T),
"mean_SOIL_ppb"=mean(SOILf_Val, na.rm=T),
"mean_SVR_km"=mean(SVR_Val, na.rm=T))%>%
as.data.frame()%>%
separate(SiteName, into = c("SiteName", "year"), sep="_")%>%
separate(year, into = c("year", "season"), sep=" ")
head(mean_by_season)
## SiteName year season mean_ammNO3_ppb mean_ammSO4_ppb mean_EC_ppb
## 1 Agua Tibia 2011 autumn 0.7042493 1.236351 0.2399840
## 2 Agua Tibia 2011 spring 1.3284576 1.929723 0.1828920
## 3 Agua Tibia 2011 summer 0.9324993 2.694624 0.2530929
## 4 Agua Tibia 2011 winter 0.6691160 0.510433 0.1339632
## 5 Agua Tibia 2012 autumn 0.6549057 1.049741 0.1992536
## 6 Agua Tibia 2012 spring 2.0285125 2.128372 0.2105500
## mean_OMC_ppb mean_SOIL_ppb mean_SVR_km
## 1 1.5392304 0.4572863 137.75418
## 2 1.7740872 0.8754060 93.61895
## 3 2.2604336 0.5468550 84.92630
## 4 0.8716737 0.2251410 191.02952
## 5 1.2527293 0.4292511 141.15836
## 6 1.7745840 0.7293125 82.39903
Creating Plots Here is our reference graph:
mean_by_season%>%
filter(SiteName=="Fresno")%>%
ggplot(aes(year, mean_SVR_km, fill=season))+
geom_col(position = "dodge")
Here is the shiny app:
ui <- dashboardPage(
dashboardHeader(title="Seasonal Metrics"),
dashboardSidebar(disable = T),
dashboardBody(
tags$head(
tags$style(HTML(".box{border-radius: 20px; box-shadow: 5px 5px 10px #888888;}"))),
fluidRow(
class="text center",
box(
selectInput("site", "Select Site Name", choices = unique(visibility_tidy_mdy$SiteName))
),
box(
selectInput("y", "Visibility Metric",
choices = c("Standard Visual Range"="mean_SVR_km", "Ammonium Sulfate" = "mean_ammSO4_ppb",
"Ammonium Nitrate" = "mean_ammNO3_ppb", "Elemental Cabon" = "mean_EC_ppb",
"Organic Mass Cabon" = "mean_OMC_ppb", "Dust" = "mean_SOIL_ppb"),
selected = "mean_SVR_km")
),
box(
plotOutput("plot", width = "1400px", height = "400px"), width = 12
)
)
)
)
server <- function(input, output, session) {
session$onSessionEnded(stopApp)
output$plot <- renderPlot({
mean_by_season%>%
filter(SiteName == input$site)%>%
ggplot(aes_string("year", input$y, fill="season"))+
scale_fill_manual(values=c("#00CCE6", "#9BEB01", "#E6E500", "#EB7500"))+
geom_col(position = "dodge", color="black")+
theme(axis.title = element_text(size = 16),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.background = element_rect(fill = "#F1F1F1"))+
labs(x="Year",
y=paste(str_split_1(input$y, "_")[2], str_split_1(input$y, "_")[3], sep=" "),
fill="Season")
})
}
shinyApp(ui, server, options = list(launch.browser = TRUE))